################################################################################################

library(parallel)
library(data.table)
library(Seurat)
library(optparse)

## 很关键
Sys.setenv(RETICULATE_PYTHON = "~/miniconda3/envs/magic_v1/bin/python3")
library(reticulate)
library(Rmagic)
# repl_python()

##########################################################################################

option_list <- list(
    make_option(c("--clustern"), type = "character"),
    make_option(c("--input_dir"), type = "character"),
    make_option(c("--output_dir"), type = "character"),
    make_option(c("--rna_data_file"), type = "character"),
    make_option(c("--qc_cell_file"), type = "character"),
    make_option(c("--atac_data_file"), type = "character")
)

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

clustern <- gsub( "cluster" , "" , opt$clustern)
input_dir <- opt$input_dir
rna_data_file <- opt$rna_data_file
qc_cell_file <- opt$qc_cell_file
atac_data_file <- opt$atac_data_file
output_dir <- opt$output_dir

###########################################################################################

if(1!=1){

	clustern <- 0
	input_dir <- "~/20231121_singleMuti/results/cluster_all_result/cluster0"
	output_dir <- "~/20231121_singleMuti/results/cluster_all_result/cluster0/magic"
	rna_data_file <- "~/20231121_singleMuti/input/testis_combined.Rdata"
	qc_cell_file <- "~/20231121_singleMuti/results/cluster_all_result/cell_fragment/testis_merge_all_qc_barcode_new.txt"
	atac_data_file <- "~/20231121_singleMuti/results/cluster_all_result/cluster0/count/cluster0_testis_merge_peak_barcode_matrix.Rdata"

}

dir.create(output_dir , recursive = T)

################################################################################################
## 所有质控合格的细胞
cell_all <- data.frame(fread(qc_cell_file,header=T))

## single cell RNA data ##
a <- load(rna_data_file)  

## counts数据，三批合并后的
b <- load(atac_data_file)

################################################################################################
## 替换atac的细胞名
cell_all$V1 <- paste0( sapply(strsplit(cell_all$cellID , "_" , 1) , "[" , 1) , "_" , cell_all$V1 )
rownames(cell_all) <- cell_all$cellID
colnames(cluster_testis_count) <- cell_all[colnames(cluster_testis_count),]$V1

################################################################################################
## 读入atac的counts矩阵
cluster_count_seu <- CreateSeuratObject(cluster_testis_count, project = "CreateSeuratObject")
cluster_count_seu <- NormalizeData(object = cluster_count_seu)
cluster_count_seu <- ScaleData(object = cluster_count_seu)

## 对counts矩阵进行填补
cluster_count_magic <- magic(cluster_count_seu, knn=10, t=6, npca=20)

## 释放内存
rm(cluster_testis_count)
#save(cluster_count_magic,file=paste0(input_dir,"/cluster",clustern,"_peak_magic.Rdata",sep=""))

################################################################################################
## 转录组测到的基因
scrna_gene <- rownames(testis_combined@assays$RNA@counts)
#scrna_gene <- read.table(paste0(input_dir,"/scRNA_gene.txt",sep=""),header=FALSE)  ## single cell RNA genes ##
cluster_overlap_effc <- fread(paste0(input_dir,"/annotation/cluster",clustern,"_peak_annoatation.txt",sep=""),header=FALSE)
cluster_overlap_effc <- subset( cluster_overlap_effc , V14>=250 & V11 %in% scrna_gene )
cluster_overlap_effc$V4 <- gsub("_","-",cluster_overlap_effc$V4)

## 提取质控合格的细胞
cell_in_rnam <- data.frame(barcode=colnames(testis_combined@assays$MAGIC_RNA))
cell_in_rnam <- merge(cell_in_rnam,cell_all,by.x="barcode",by.y="V1")
rna_data <- testis_combined@assays$MAGIC_RNA[,cell_in_rnam$barcode]
colnames(rna_data) <- cell_in_rnam$barcode

atac_data <- cluster_count_magic@assays$MAGIC_RNA
rna_data <- rna_data[,colnames(atac_data)]

save(rna_data,atac_data,cluster_overlap_effc,file=paste0(output_dir,"/cluster",clustern,"_magic_for_correlation.Rdata",sep=""))


